A new shear creep damage model for rock masses after considering initial damage

For a long time, one of the important safety problems in open-pit mines is the stability of a large number of high slopes with gently inclined soft interlayer. Rock masses formed after long geological processes generally have some initial damage. Mining works also cause varying degrees of disturbance and damage to rock masses in the mining area during the mining process. This phenomenon means that accurate characterization of the time-dependent creep damage for rock masses under shear load is necessary. The damage variable D is defined based on the spatial and temporal evolution laws of shear modulus and initial level of damage for the rock mass. In addition, a coupling damage equation between the initial damage of the rock mass and shear creep damage is established based on Lemaitre’s strain equivalence assumption. Kachanov’s damage theory is also incorporated to describe the entire process of time-dependent creep damage evolution for rock masses. A creep damage constitutive model that can reasonably reflect the actual mechanical properties of rock masses under multi-stage shear creep loading conditions is established. This takes into account multi-stage shear creep loading conditions, instantaneous creep damage during the shear load phase, staged creep damage and factors influencing the initial damage of rock masses. The reasonableness, reliability and applicability of this model are verified by comparing the results of the multi-stage shear creep test with calculated values from the proposed model. As opposed to the traditional creep damage model, the shear creep model established in this present study takes into account the initial damage of rock masses and can describe the multi-stage shear creep damage characteristics of rock masses more convincingly.


Introduction
Time-dependent creep behavior is an important indicator in the evolution of mechanical properties for rocks and greatly affects the long-term stability of rock masses [1][2][3][4]. Slipping and collapsing slopes, mined-out areas and dam foundations are usually caused by rheological behavior [5,6]. Therefore, it is of great significance to study the creep properties of rock masses to reveal their creep mechanisms and damage criteria so as to facilitate the evaluation of longterm stability in rock engineering. Since Kachanov put forward the damage theory in 1985 (the exponential function used to describe the time-dependent damage of rock mass can better characterize the time-dependent damage of rock mass) [7-9], creep properties of rocks have been thoroughly studied over the last 70 years, with a large number of mechanical tests [10][11][12][13][14] carried out and constitutive models [15,16] established. Various constitutive creep models such as Maxwell, Burgers, Bingham and Nishihara have been used to simulate the mechanical behavior of rock masses under long-term loading. Major improvements have been made for these traditional creep models to improve their accuracy and applicability [17,18].
Element models based on mechanical elements have undergone extensive development in the field of rock mechanics research [19]. Wang et al. [1] established a three-dimensional creep damage model based on the traditional Burgers creep model by considering creep damage with respect to time and the effect of disturbance damage in the spatial dimension. Liu et al. [20] introduced damage variables in the Nishihara model to characterize the damage process of viscosity coefficients and studied the nonlinear creep properties of soft rocks. Zhou et al. [21] proposed a creep damage constitutive model based on time-fractional derivatives by replacing Newtonian dampers in the classical Nishihara model. Zhang et al. [22] proposed a Nonlinear Burgers Strain Softening (NBSS) model and a numerical calculation method to describe the plastic zone and deformation law for the surrounding rock. Di et al. [23] put forward a fractional viscous-elastic-plastic model, showing that both creep and recovery curves follow a power-law model. Wang and Cai [24] modified the viscosity coefficients of Maxwell components and proposed a time-to-failure (TtoF) creep model. Xu et al. [25] proposed a novel creep model connecting Maxwell bodies, Kelvin bodies and nonlinear visco-plastic bodies, which described the entire creep process and transverse isotropic characteristics of phyllite. However, the mathematical complexity and redundant parameters of these multi-component systems make practical applications difficult.
Fractional calculus has been applied to describe the viscoelastic and Visco-Plastic deformation behavior of materials [26][27][28][29]. Welch et al. [30] proposed a fractional calculus-based constitutive model to describe the elastic viscosity time-dependent deformation of polymeric materials. Tang et al. [31] proposed a nonlinear creep model based on fractional derivatives and damage mechanics theory to describe rock characteristics at different creep stages. Chen et al. [32] proposed a time-dependent damage constitutive model based on fractional calculus and damage mechanics theory. Xu et al. [33] investigated a fractional constitutive relationship model represented by a fractional element network. This showed the effectiveness of the Interior Point Method (IPM) in parameter inversion of the fractional viscoelastic model. Liu et al. [34] proposed a variable-order fractional-order derivative model with damage taken into consideration. The visco-plastic response for parallel combinations of variable-order fractional components and plastic slides was also described. Kang et al. [35] developed a fractional nonlinear creep model using Scott-Blair fractional elements to describe the viscoelastic-plastic creep behavior of coal and characterized the accelerated creep process in terms of damage factors. The constitutive model based on fractional calculus has some limitations. It is difficult to reasonably explain the initial damage of rock mass from a theoretical perspective.
The study of creep damage in rocks has yielded a wealth of results. However, little research has been reported on the rational representation of damage variables. It is necessary to consider the initial damage factors of rocks and the rational representation of the time-dependent damage evolution for rock shear creep based on coupling damage (initial damage and shear creep damage). To address this problem, a new shear creep damage model is developed in this paper based on Kachanov's damage law and Lemaitre's strain equivalence principle.
Computational results obtained by parameter inversion are compared with experimental results to verify the feasibility and applicability of the proposed model.

Evolutionary process of deformation modulus for rock masses under multi-stage shear creep loading and definition of damage variables
Rocks naturally form pores, fissures, cracks, dislocations and other defects at different stages of time during their long geological history, thereby leading to uneven initial damage in rocks [36][37][38][39]. Current research on creep damage theory in rocks generally assumes that the rock is in an undamaged state at the initial creep stage (t = 0) [40][41][42].
Considering the initial damage of the rock mass and the definition of damage variables, the damage variable is D = D a = 1 -G b /G 1 when the rock is not damaged by shear load. The initial shear modulus of the rock shear surface is G b The shear modulus of rock mass as a whole is G 1 . Based on shear damage characteristics of the rock mass, when no initial shear load is applied to the rock, we have G = G b , D = D a ; when the rock is sheared along the shear damage surface, we have D = 1. At this time, the rock is completely fractured along the shear damage surface and the shear damage surface as a whole is represented by G = 0 without residual shear modulus. Since the initial damage to the shear surface of the rock mass does not change with shear creep load, the shear modulus of the rock mass is considered equivalent to the initial shear modulus of the shear surface. Given the initial damage to the shear surface of the rock mass and the evolutionary characteristics of the shear modulus, the damage variables of the shear surface for the rock mass under multi-stage shear creep loading are defined as: Where: G is the time varying shear modulus (secant modulus) of rock masses under multistage shear creep load.

Expressing the level of multi-stage creep damage of rock masses after considering initial damage
Based on the multi-stage creep test data of rock masses under shear load obtained from the test, the full stress-strain curve of rock masses is obtained through processing, and the shear modulus (secant modulus) under stress and time conditions is calculated. The calculation method of the time varying shear modulus value is similar to that of elastic modulus, which is the slope of secant, as shown in Fig 1. According to the time-dependent nature of multi-stage creep, the secant modulus can be divided into two parts. These are the instantaneous secant modulus G is obtained from the instantaneous shear deformation calculation and the secant modulus G iT after multi-stage creep obtained from the creep deformation calculation.
where τ is , ε is , ε iT are the shear stress, instantaneous strain and creep end strain of the rock mass when subject to loading in stage i.
In order to further improve the reliability of the time-dependent damage calculation in Eq (1), the initial damage was considered when defining damage variables. However, the coupling effect of initial damage on time-dependent damage was not considered. According to literature, most of the current studies on multi-stage shear creep models for rock masses are conducted based on the element model. This only considers the time accumulation effect and the microscopic damage of rock masses caused by shear stress [43,44]. Multi-stage shear creep damage models that consider the effects of initial damage such as rock mass shear surface joints and fractures have not been found. However, in practice, the stability of rockwork is controlled by macroscopic defects such as joints (Fig 2) and destabilization damage is usually caused by such initial damage [45,46]. Therefore, when applying the model in practice, it is necessary to consider both initial damage caused by natural or anthropogenic factors and microscopic time-dependent damage due to rheology.

PLOS ONE
To obtain the coupling damage equation for time-dependent creep damage under the influence of the initial damage at the shear surface of the rock mass, we can refer to the elastic-plastic mechanics of the rock mass. Several parameters can be obtained, such as the timedependent relationship under shear load between the effective shear modulus at the shear surface of the rock mass, the initial effective shear modulus at the shear surface of the rock mass, the time-dependent shear modulus(secant modulus) during the destruction of the rock mass along the shear damage surface and the initial effective shear modulus of the rock mass as a whole. These can be written as: Based on Lemaitre's strain equivalence assumption and taking into account the initial damage to the shear surface, we can obtain: Substituting Eq (6) from above into Eq (5), we can obtain the coupling damage equation for time-dependent creep damage that has taken into account the initial damage of the rock mass.
From the above equation, it can be seen that the initial damage of the rock shear surface is D a = 0 with coupling damage represented by D c = D at . The time-dependent creep damage of the rock shear surface under shear load is D at = 0 with coupling damage being D c = D a . These results show that the coupling damage expression is applicable to the study of mechanical properties for shear creep damage coupling under the influence of initial damage.
Based on Lemaitre's strain equivalence assumptions, the coupling damage expression obtained in [47] can be written as: Based on the formula proposed in this paper to calculate the coupling damage of the rock mass shear surface and the coupling damage expression for the rock mass obtained in [47], the coupling damage variables calculated are shown in Fig 3. It is not difficult to conclude that the coupling damage variables obtained in this paper are the smallest. As the number of macroscopic and microscopic fractures in the rock mass increases during the creep process, initial damage accounts for a reduced proportion of total coupling damage. When creep damage reaches 0.6, the proportion of the initial damage (0.05) in the coupling damage is only 0.0082, which is consistent with the shear creep damage law for the rock mass. The coupling damage variable in [47] is significantly larger than the coupling damage variable in this paper. Under the condition of constant initial damage for the rock mass, the rate of change for coupling damage Ḋ c as creep damage increases is always constant. The initial damage is not reduced when the creep damage increases. This shows that the coupling damage method proposed in [47] is not applicable to the expression for shear creep coupling damage (initial damage and creep time damage) in rock masses.
According to the theory of rock damage mechanics, the expression for initial damage can be written as: From the above equation, when the rock mass is sheared along the shear surface, the shear modulus of the shear surface is G b = 0 and shear slip of the rock mass occurs along the shear damage surface.

PLOS ONE
Substituting Eq (9) from above into Eq (7) and combining this with Eq (1) yields the coupling damage equation for initial damage and time-dependent creep damage in the shear surface of the rock mass.

Constitutive equations for multi-stage shear creep damage based on Kachanov's damage law and Lemaitre's strain equivalence principle
Based on Lemaitre's strain equivalence assumptions, the corresponding shear damage constitutive equation is given by: wheret is the effective shear stress at the shear surface andG is the effective shear modulus at the shear surface.
According to Lemaitre's strain equivalence principle, the relationship between the shear surface strain and the time-dependent strain for shear load at the i th -stage is: The above equation shows that when D it (t) = 0, the shear surface has no shear creep damage and the effective modulus isG ¼ G is . When D it (t) = 1, the effective shear modulus of the shear damage surface isG ¼ 0, indicating that Eq (12) is consistent with the general law of shear creep damage. However, for multi-stage shear creep loading conditions, it is not only the instantaneous creep damage that should be taken into consideration during the shear stress loading stage. The stage creep damage and the initial damage of the rock mass should be considered as well. Therefore, in order to establish a constitutive equation that can reasonably reflect the actual mechanical properties of a rock at its shear surface, the above equation is modified in this paper and the modified expression for effective shear modulus is: It can be seen from the above modified effective shear modulus equation that when D it (t) = 0, there is no shear stress loading and no load damage occurs on the shear surface, with the effective shear modulus calculated asG ¼ G b . When D it (t) = 1, the effective shear modulus can be calculated asG ¼ G iT � 0 which is in line with actual working conditions under the action of a multi-stage shear creep load. Substituting Eq (13) into the strain expression in Eq (12), the modified shear damage constitutive equation of the rock mass is obtained as follows: where ε is is the strain after shear loading at the i th -stage. Based on the damage variables proposed by Kachanov, the damage evolution equation is [48]: where δ and ξ are material constants determined from the creep test. When integrating Eq (15) with the initial conditions when the rock mass is completely destroyed along the shear surface, it is found that D e = 1. In this case, we can get the expression for the end time t u of shear creep behavior for the rock mass, which is written as: Combining Eqs (15) and (16), the evolution equation of Kachanov's creep damage law can be obtained as follows: Under the action of the last creep load, when the shear plane undergoes accelerated creep failure, the rock mass enters the accelerated creep damage stage. The parameters ξ s of the Kachanov damage law have certain limitations in modifying the creep characteristics of rock masses in the accelerated creep damage stage. To more accurately describe the damage characteristics of rock masses in the accelerated creep stage, the Kachanov damage law correction coefficient η is introduced, and a piecewise function is constructed to reasonably express the time varying shear modulus damage of the shear surface.
where t u is the critical time when creeping ends for the shear surface and the rock is just about to enter the accelerated creep phase; ε s and G s are the shear strain and shear modulus at the beginning of the accelerated creep phase on the shear surface respectively; ξ S is the Kachanov damage law parameter for the accelerated creep phase; t a ð Þ ¼  (14), we can obtain the shear surface damage constitutive equation of the rock mass under loading at the i th -stage as follows: As can be seen from Eq (18), when ε < ε s , z is an important parameter reflecting the creep damage trend of the rock shear surface; when ε > ε s , z s and η are important parameters reflecting the creep damage trend of the rock shear surface at the accelerated creep damage stage. A better understanding of the influence for parameters z, η and z s on the creep damage characteristics of the rock shear surface at the accelerated creep damage stage is required. As such, the sensitivity of the parameters z, η and z s to the shear creep damage characteristics of the rock mass is analyzed. When ε < ε s , to better reflect the influence of initial damage on the time-dependent creep damage, the parameters are taken to be the initial creep moment. In order to better reflect the influence of initial damage on aging creep damage, the parameters are recorded based on the initial creep moment, i.e. G is = G b ; when ε > ε s , the variables D a = 0.06, G b = 0.296Gpa, G is = 0.234Gpa, G e = 0.102, G s = 0.186, t u = 18 and t f = 24 in Eq (18) are chosen. The evolution of the damage variables D it (t) with creep time was calculated for different parameters z, η and z s .
As shown in Fig 4, when t = 0, the damage variables are represented by the initial damage of the rock mass. When ξ and ξ s > 0, as creep time increases, the time-dependent damage evolution rate Ḋ it of the rock mass keeps increasing. The larger the values of z and z s , the more significant the growth trend is. When ξ and ξ s < 0, as creep time increases, the time-dependent damage evolution rate Ḋ it of the rock mass keeps decaying. The smaller the value of z and z s , the more significant the decay trend is. This indicates that both creep damage parameters and time-dependent damage variables of the rock mass show a certain sensitivity. It may be the case that changes in the creep damage parameters cannot satisfy the evolution trend of the accelerated creep phase for the rock mass. For this, an accurate representation of the accelerated creep phase can be achieved by adjusting the introduced Kachanov damage law correction coefficients, thereby making it more widely applicable. This indicates that both the creep damage parameters z and z s show a certain sensitivity to the time-dependent damage variable D it (t) of the rock mass. It may also be the case that changes in z s cannot satisfy the evolution trend of the rock mass in the accelerated creep phase ε-t. To resolve this, an accurate representation of the accelerated creep phase can be achieved by adjusting the introduced Kachanov damage law correction coefficient η, thereby making it more widely applicable.

Model validation
In order to verify the broad applicability of the proposed creep model, this paper uses experimentally obtained multi-stage shear creep test data of the rock mass under the influence of impact disturbance factors.

Shear creep test
A self-developed WJSZ-100 microcomputer servo controlled rainfall seepage and blasting vibration coupling simulation shear rheology test system of soft rock was adopted in the test. The test system consists of four modules: compression shear subsystem, seepage subsystem, excitation subsystem, and data monitoring and acquisition subsystem. Some parameters of the test system are as follows: Axial load control range: 0~100kN; Shear load control range: 0~100kN; Normal loading rate: 0.001~10mm/min; Shear load loading rate: 0.001~10mm/min; Load holding time: 2000h; Range of normal displacement sensor: 0~30mm; Range of tangential displacement sensor: 0~30mm; Measurement accuracy of normal/tangential displacement sensor: 0.001mm.
Description of shear creep test principle: The horizontal shear load loading device on the right side of the frame is used to apply shear load to the sample (the loading direction is from right to left), and the reaction column is set on the opposite side of the horizontal static load loading device (as shown in Fig 5). The direction of the normal stress is axial (vertical direction). The normal stress loading device is set in the upper part of the frame, the loading head is set as a ball and the upper part of the shear box is located the same centerline. And the normal stress loading direction is vertical (from top to bottom). The loading of the normal stress in the shear creep test plays a role in fixing the shear box and preventing the shear box from deviating along the normal direction during the loading of the shear creep load. The displacement

PLOS ONE
sensor and strain gauge of the test system were used to record the strain deformation of the sample every 0.5s. After the test, the shear force and normal pressure were removed successively, take out the sample and process the test data.
The tested rock samples were taken from the mud shale of Jinding Mine in Emeishan, Sichuan Province (China), and the test samples were taken from the same shale. Samples were numbered A 1 , A 2 , A 3 and A 4 . A multi-stage shear creep test was carried out on the samples (Fig 5). During the test, the axial pressure was kept at 1 MPa and the shear loading process was divided into six stages: 0.81 MPa, 1.21 MPa, 1.61 MPa, 2.02 MPa, 2.42 MPa and 2.82 MPa. Each stage lasted for 24 hours and the impact disturbance load was applied at the 17 th hour of the creep stage. The impact energies of samples A 1 , A 2 , A 3 and A 4 were 0 J, 1.34 J, 2.36 J and 3.67 J respectively. And process specification of samples testing. Related test parameters for the mud shale are shown in Table 1. The multi-stage shear creep curve was obtained using multistage shear creep test data (Fig 6).
It can be seen from Fig 6 that under the condition of multi-stage shear creep load, with the increase of shear creep load, the deformation of rock mass gradually increases, and there is an obvious accelerated creep stage before the complete failure of rock mass along the shear plane.

PLOS ONE
In the process of decelerating creep stage and stable creep stage, the shear surface is continuously damaged to different degrees. With the increase of cumulative damage on the shear plane, the shear failure occurs within a few minutes after the specimen enters the accelerated creep stage. And the creep rate in the steady creep stage is significantly higher than the steady creep rate at low shear stress. The impact disturbance energy has a certain influence on the creep characteristics of shale. Under the same normal stress, with the increase of impact disturbance energy, the creep deformation rate of the sample increases gradually, and the ductile failure characteristics increase. When the impact disturbance energy is small (1.34 J), the damage to the specimen is small, and the specimen is damaged under the sixth order shear stress; When the impact disturbance energy is large (2.36 J, 3.67 J), the damage degree of the sample

PLOS ONE
will be increased, and the sample will be damaged under the fifth level shear stress. With the increase of the impact disturbance energy, the creep failure time of the rock mass under the fifth level shear stress will be shortened. The influence of shear creep load and impact disturbance on the sample is mainly the change of physical and mechanical properties of the sample, which promotes the destruction of the cemented filler between mineral particles in the shale, the rolling and shearing of mineral particles, the initiation, and expansion of cracks in the sample, and thus increases the damage degree of the sample; With the increasing of shear creep load and impact disturbance energy, the cracks inside the shale sample are constantly damaged and deteriorated, the strength of the sample decreases. Finally, the sample undergoes accelerated creep deformation under the last shear load.

Creep damage model validation
According to the data in Table 2, the test results in Fig 6 are curve-fitted using Eq (19) and the fitting results are shown in Fig 7. It can be seen that the calculated values for this model are in good agreement with the experimental data, R 2 > 96%. The results show that the creep damage model proposed in this paper is reasonable and reliable under multi-stage shear creep load conditions and can effectively reflect the deformation and failure characteristics of rock masses. Table 3 shows inversion results for sample parameters ξ , ξ s , η.

Discussion
Based on the Lemaitre strain equivalence hypothesis and Kachanov damage theory, a new constitutive model of rock mass shear creep damage has been established in this paper, and its Because the model proposed in this paper is based on rock mechanics and related damage mechanics theory, which has certain theoretical rigor and feasibility. And the secant modulus method is used to define the damage variable to avoid the difference of related mechanical parameters of different kinds of rocks. Moreover, the relevant parameters (z and z s ) of the model are within the exponential term of the improved Kachanov damage law exponential function. And the Kachanov damage law correction coefficient η is introduced, which can further modify the function curve in the accelerated creep stage. So it has a strong regulatory effect on the rock shear creep ε-t function curve. The uncertainty corresponding to its applicability lies in the adjustment limit of the model parameters (z 、η and z s ) in the shear creep function curve of rock masses, which needs further study. The establishment of the shear creeps model considering the initial damage of rock masses can provide a theoretical basis for the establishment of an early warning system for the actual mine high slope rock landslide (sliding along the weak structural plane-weak interlayer). With the continuous deepening of the research on the creep damage model of rock masses, the future research focus is to reasonably select and improve the existing creep model of rock mass based on the specific rock mass structural characteristics of the actual mine and the relevant data of the measured regional rock masses deformation, predict the deformation rules of the actual mine rock masses, to achieve the purpose of disaster early warning.

Conclusion
The new shear creep damage constitutive model based on the spatio-temporal evolution characteristics of shear modulus, Lemaitre strain equivalent hypothesis and Kachanov damage theory has simple structure and fewer parameters compared with the element creep constitutive model, and can be well used to reveal the evolution law of multi-stage shear creep of rock mass. The main conclusions are as follows: 1. A coupling damage expression for initial damage of the rock mass and shear creep damage was obtained based on Lemaitre's strain equivalence assumptions. By relying on the representation of initial damage for the rock mass, the reasonableness and applicability of the proposed coupling time-dependent damage equations were analyzed using a comparative analysis method based on Kachanov's damage theory.
2. Based on the actual working conditions for multi-stage shear creep damage, a novel and reasonable shear creep damage constitutive model was established. It takes into consideration instantaneous creep damage during the shear loading stage, the stage creep damage and the initial damage factors of the rock mass.
3. Under the condition of constant initial damage for rock masses, the rate of change for coupling damage Ḋ it as creep damage increases is always constant in [47]. The proportion of initial damage is not reduced when the creep damage increases in [47]. This shows that the coupling damage method proposed in [47] is not applicable to the expression for shear creep coupling damage (initial damage and creep time damage) in rock masses. The feasibility and applicability of the proposed shear creep model were demonstrated based on data related to multi-stage shear creep that were obtained via testing.
4. The research results of this paper show that the spatial and temporal evolution of shear modulus can help to realize a quantitative damage expression under the rock mass shear creep law. This applies even for the case of different disturbance factors and facilitates the study of the influence of different disturbance factors on multi-stage shear creep damage of rock masses from the perspective of damage quantification.